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We propose an efficient iterative method to solve the mixed Stokes-Dracy model 
for coupling fluid and porous media flow. The weak formulation of this problem 
leads to a coupled, indefinite, ill-conditioned and symmetric linear system of equa- 
tions. We apply a decoupled preconditioning technique requiring only good solvers 
for the local mixed-Darcy and Stokes subproblems. We prove that the method is 
asymptotically optimal and confirm, with numerical experiments, that the perfor- 
mance of the preconditioners does not deteriorate on arbitrarily fine meshes. 
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The Stokes-Darcy problem describes filtration processes that find many important ap- 
plications in porous media problems. Usually, a surface free flow of a liquid is modeled 
by Stokes equations and the flow confined in the porous media is governed by Darcy 
equations. The interaction of the local models is commonly handled through the Beavers- 
Joseph-Saffman (BJS) interface conditions, cf. [271 |2"T] . 

Recently, there has been active research on the mathematical and numerical analysis 
for this model. The weak formulation of this problem is generally obtained by coupling 
the usual velocity-pressure mixed formulation in the Stokes domain with either the primal 
formulation (i/ 1 -approach, [15J) or the mixed formulation (H(div)-approach, [22l [T8] ) in 
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the Darcy domain. In both cases, optimal iterative methods are of crucial importance to 
solve efficiently the discrete Stokes-Darcy model since the corresponding linear systems 
of algebraic equations are indefinite and ill-conditioned. Several optimal iterative solvers 
such as the Dirichlet-Neumann or Robin-Robin domain decomposition methods PUJ dU 
[T2"l [T3J [TJ], multi-grid methods (SJ EJ 125] have been proposed for the model based on 
an if ^approach in the Darcy domain. The common feature for most of these methods 
consists in decoupling the global model in such a way that, only independent Stokes and 
Darcy subproblems are involved in the iterative process. To the authors' knowledge, the 
only solution procedure for the H(div)-conforming Darcy flux approach was proposed in 
[2]. In this paper, the problem is written as a global saddle point problem and a solver 
is implemented for the scheme using an inexact Uzawa technique relying on an expensive 
preconditioner. Our purpose here is to devise a decoupled preconditioning strategy that 
allows to apply existing optimized solvers to each local model independently. 

Recently, a systematic way to obtain convergent finite element schemes for the Darcy- 
Stokes flow problem by combining well-known mixed finite elements that are separately 
convergent for the H(div)-Darcy formulation and the Stokes problem was proposed in 
|24j . We take advantage here of a fluid-to-pressure (FtP) operator to reinterpret in this 
formulation the Darcy system as a nonlocal boundary condition for the Stokes problem. 
The corresponding discrete equations are written in terms of a symmetric and indefi- 
nite linear operator that enjoys the same spectral properties of the local discrete Stokes 
problem. 

Many different iterative methods for solving the saddle point problems that result from 
the finite element discretization of the Stokes equations are known. There are, for instance, 
many variants of the so-called inexact Uzawa methods. Block triangular preconditioners 
for saddle point operators has also been discussed by many authors. An example of such 
preconditioners has been introduced in j5] by Bramble and Pasciak. In this strategy the 
original saddle point system is premultiplied by a block triangular operator and then, 
the resulting positive definite system is solved by a preconditioned conjugate gradient 
method. A possible difficulty is that this approach requires a proper choice of a critical 
scaling parameter to obtain a positive definite operator. Finally, we mention the block 
diagonal preconditioners for the minimal residual method (MINRES), cf. [T6] and 
the references therein. A comparative study of representing methods from each of these 
three classes is considered in [26]. The inexact Uzawa method is not feasible in our case 
because of the nonlocal character of the discrete flux-to-pressure operator Cg appearing 
in the principal block of our saddle point problem ( 14. ip . The conclusion in [2E] is that the 
preconditioned MINRES method may be slower than the Bramble-Pasciak method but it 
is more robust (it even converges without preconditioning) and it is parameter free. 

Applying a preconditioned MINRES method in our case requires, at each iteration 
step, the solution of two local problems: a vector Laplace equation in the fluid and the 
mixed formulation of the Darcy problem in the porous media. This saddle point problem 
in the Darcy domain is again solved with a preconditioned MINRES method. The global 
algorithm has then the structure of an outer-inner MINRES iteration process. Thus, for 
this method a stopping criterion (tolerance parameter) for the inner iteration is needed. 
The preconditioner for the inner MINRES may be constructed by using techniques from 
[31 ED] • We use here the nodal auxiliary space preconditioning technique introduced in [2D] 
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by Hiptmair and Xu to solve H(div)-elliptic problems. With this choice, our decoupled 
iterative process consists in two nested MINRES methods whose preconditioners only 
require the solution of several second-order if ^elliptic problems in the Stokes and the 
Darcy domains. Standard multigrid techniques or domain decomposition methods can 
then be applied to reduce the computational effort. Theoretical analysis and numerical 
experiments show the optimality and efficiency of the proposed decoupled iterative solver. 

The rest of the paper is organized as follows. In Section [2] we summarize the results 
of [23] introducing the model problem, the variational formulation, general conditions 
for convergence of a Galerkin discretization and examples of spaces leading to convergent 
methods. A reinterpretation of the continuous formulation in terms of a Fluid-to-Pressure 
operator and the derivation of its discrete counterpart are presented in Section [3] We take 
advantage of the equivalent formulation of the discrete Stokes-Darcy problem to deduce, 
in Section HI a decoupled iterative solver based on a preconditioned MINRES method. 
Finally, numerical experiments are reported in Section |5j 

Notation and background. Boldface fonts will be used to denote vectors and vector 
valued functions. Also, if if is a vector space of scalar functions, H will denote the space 
of IR d valued functions whose components are in H, endowed with the product norm. 

Given an integer m > 1 and a bounded Lipschitz domain O C M. d , (d = 2,3), we 
denote by || ■ || ifm(o) the norm in the usual Sobolev space H m (0), cf. [T]. For economy 
of notation, (•, -) stands for the inner product in L 2 {0) and || • \\o is the corresponding 
norm. We recall that H l / 2 (dO) represents the image of H l (0) by the trace operator. 
Its dual with respect to the pivot space L 2 (dO) is denoted H~ 1 ' 2 (dO). For definition 
and basic properties of the spaces H(div, O) and H(curl, O), we refer to [IS]. We will 
denote by H (div, O) the subspace of fields from H(div, O) with zero normal trace on 
the boundary dO and by H (curl, O) the subspace of fields from H(curl, O) with zero 
tangential trace on dO. 

For k > 0, Pfc(O) will denote the space of d— variate polynomials of degree not greater 
than k defined on the set O C M. d with non-trivial interior. Finally, at the discrete level, 
the letter h (with or without geometric meaning) will be used to denote discretization. 
The expression a <b will be used to mean that there exists C > independent of h such 
that a < Cb for all h. Similarly, we write a ~ b when there exist constants C > c > 
independent of h such that ca < b < Ca. 

Let us consider a linear operator A h : — > X^ acting between a finite dimensional 
subspace X^ of a Hilbert space X C L 2 (0) and its dual X^. Assume that we have chosen 
a basis in Xh and that the coefficients of Uh, Vh € Xh in this basis are given by u, v G M n , 
where n is the dimension of Xh. We define the matrix realization A h e R nxn of A h by 

{A h u, v) 2 = {A h u h , v h ) x * h xx h Vw h , v h e X h , (1.1) 

where (•, stands for the Euclidean scalar product in MJ 1 . Moreover, if I h : X^ — > X' h is 
the Riesz operator given by 

(I h u h ,v h ) x * h xx h = (u h ,v h ) Vu h ,v h G X h , 
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then, the corresponding matrix realization M e M nxn , usually referred to as the mass 
matrix, is defined by 

(M h u, v) 2 = (I h u h , v h ) x *xx h Vm a , v h G X h . 

2 Statement of the problem and discretization 

Let us consider a domain fi C M. d (d = 2 or d = 3) with polyhedral Lipschitz boundary. 
We assume that fi is subdivided into two subdomains by a Lipschitz polyhedral interface 
E. The subdomains are denoted Slg and fin (S stands for Stokes and D for Darcy). We 
also denote Tg := <9fig \ E and := dQj) \ E. The normal vector field n on dQ is chosen 
to point outwards. We also denote by n the normal vector on E that points from Q s to 













E^^ — 















2.1 Variational formulation 

In the region fig, the fluid flow is assumed to satisfy the Stokes system 

-div(2i/e(u s ) -psl) = fs, divu s = in fig, (2.1) 

where I is the identity in M. d and e(u s ) := |(Vug + V T u s ) is the deformation tensor, 
v > is the kinematic viscosity and fs is the external body force. In the porous region 
fio, the governing equations are given by the following Darcy system 

K- 1 u D + Vp D = 0, divu D = / D infi D , (2.2) 

where /d is the source (or sink) term and the hydraulic conductivity tensor of the porous 
medium K(x) is symmetric and uniformly bounded and positive definite, i.e., 

< hl^l 2 < K(x)£-£ < k 2 \£\ 2 fora.e. x€fi D , V^l d , 

for some constants k 2 > k\ > 0. On the outer boundaries we consider the homogeneous 
(essential) boundary conditions 

u s = on Tg, and u D ■ n = on T D , (2.3) 



4 



and on the interface between the fluid and porous media regions we impose conditions 
ensuring mass conservation, balance of normal forces and the Beavers- Joseph- Saffman 
condition [H [27], 

u s n = u D n, 2z/£(u s )n - p s n + K7r t u s = -p D n on E, (2.4) 

where 7r t w := w — (w • n)n and k is a positive and bounded function depending on K, 
z/, and an experimentally determined friction constant, cf. (H El [TO] . 

Because of the mass conservation condition across E, the homogeneous Dirichlet 
boundary condition for us on Ts and the incompressibility condition in the Stokes domain, 
we can easily show that 

/ /d = (2.5) 

is a necessary condition for existence of solution. The pressure field is defined up to an 
additive constant. We will normalize it by imposing that 

Pd = 0. 

For the velocity field, we will use the space 

X := {u = (u s , u D ) G Hg(fi s ) x H D (div, fi D ) : u s ■ n = u D ■ n on E} C H (div, f2), 
where 

nl(n s ) := {uGH 1 ^) : u = onTg} (2.6) 
H D (div,fi D ) := {u6H(div,n D ) : u-n = on T D }. (2.7) 

The space X will be endowed with the product norm. The space for the pressure field is 
Q := L 2 (n s ) x L 2 (n D ), where L 2 (O) := {p G L 2 (0) : (p, l) c = 0}. The pressure field 
is represented as p := (ps,Pd) £ Q- Adding an appropriate constant in a postprocessing 
step, the normalization condition (p, l)n D = can be modified to (p, 1)^ = 0. The space 
Q is endowed with the corresponding product norm. 

We consider four bilinear forms, two in the Stokes domain and two in the Darcy 
domain: 



as(u s ,usj 
a D (u D ,u D ) 

frs(u s ,<?s) 
6 d (u d ,<?d) 



2u(e(u s ), e(-v s ))n s + (/«7r t us, 7r 4 v s ) s , (2.8) 

(KTHi^Vd)^, (2.9) 

(divu s ,g s ) ns , (2.10) 

(divu D ,g D ) QD . (2.11) 



These bilinear forms are combined to build the diagonal bilinear form of the mixed problem 
a : X x X — > R, given by 

a(u, v) := a s (u s ,u s ) + a D (u D ,u D ), 
as well as b : X x Q — > R given by 

6(u, q) := 6 s (u s , qs) + &d(u d , ?d). (2.12) 
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A well posed variational form of the Darcy- 
find (u, p) G X x Q such that 



Stokes problem (cf. [24"l Proposition 2.3]) is: 



a(u, v) - b(v,p) = (f s , v s )n s Vv G X, 

(2.13) 

b(u,q) ={f D ,q D ) nD Vg G Q. 

2.2 The discrete problem 

We start by creating shape-regular triangulations {T^}h and {T£}h of Qs and Qj) re- 
spectively, consisting of triangles (tetrahedra in the three dimensional case) of diameter 
not larger than h. The triangulations create two inherited partitions of S, respectively 
denoted Eg and S^. Let us consider finite dimensional subspaces of piecewise polynomial) 
to approximate velocity and pressure in the Stokes domain 

V\Q S ) C H\Q 8 ), L h (Q s ) c Lg(fts), L h (Q s ) = L h (Q s ) © P (O s ), 

as well as in the Darcy domain 

H fe (o D ) c H(div, n D ), L h (n D ) c Lg(n D ), L h (n D ) = L h (n D ) © p (o D ). 

We also need to consider the spaces 

V£(fi s ) := H h (fi s ) n H s (fi s ), H£(ft D ) := H fe (fi D ) fl H D (div, fi D ), (2.14) 

as well as the discrete spaces of normal components on £, namely, 

$ s l := {u h n : u h G V§(fl s )} C L 2 (S), 
:= {u„-n : u A G H^(fi D )} C L 2 (E). 

We will assume that contains at least the space of piecewise constant functions on 
and denote by the L 2 (E)-orthogonal projection onto 

The method we are proposing is a Galerkin discretization of the variational problem 
(I2.13P using the spaces 

X h := {u, = (u§, u£) G V^(fig) x H^(fi D ) : u£ ■ n = R h D (u h s ■ n) on £}, 
Q h := L h (n 8 ) X 4(0 D ), 

that is, we look for (u^,p^) G X ft x Q h such that 

a(u h , v ft ) - b(v h ,p h ) = (fs,v A )n 8 Vv ft G X ft , 

(2.15) 

6(u A ,g h ) =(fo,q h )n D Vq h eQ h . 

Remark 2.1. iVote £na£ X h G: X unless $g C ^ n which case (12.151) becomes a conform- 
ing Galerkin approximation of ( 12.131) . We will say that the discretization is conforming 
if $g C $d (and therefore X h G Xj and non- conforming otherwise. 
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The following result is proved in [2U Proposition 3.2]. Inf-sup conditions are written 
in terms of the spaces 



V£(fts) := V h (fl s ) n Hj(n s ), H^(fi D ) := H h (n D ) n H (div, n D ), 

which arise from the application of the discretization method to problems with homoge- 
neous boundary conditions on the entire boundary of each subdomain. 

Theorem 2.1. Let us assume that there exist a linear operator Lh : $d Hp(Qd) 

a general positive constant (3, independent of h, such that: 

sup ^ V1 ^)"s >/%J ns Vg.G/^s), (2.16) 
sup f^'^ >)9|Mk Vg fc eL*(fi D ), (2.17) 

0^u h eHg(Q D ) ll U fe||H(div,n D ) 

divH h (fi D ) C L h (fi D ), (2.18) 
3v k eVj(fl s ) s.t. (v fc • n, l) s > /3 and |K|| H i(n s) < A (2.19) 

(L>h<f>h) ■ n ||L|»0fc||H(div,n D ) < /3||0||h-i/2(s) V0 ft G (2.20) 

Taen the discrete equations (I2.15P are uniquely solvable and the following error estimate 
holds: 

||u-u fc ||x+ \\p-Ph\\n < inf ||u s - u||| H i(n s ) + inf ||u D - u£ ) || H (div,n D ) + 

u§ 6 V§(fi s ) <6H^(S1 D ) 

inf \\p - q h \\ n + \{h) (||p D - -RdPd||s + ||u s ■ n - R^(u s • n)|| s ) • 

Here X(h) = if ^ C <3>d ana 1 A(/i) < /i 1/2 otherwise. 

Let us briefly discuss the five hypotheses in Theorem 12.11 The inf-sup condition 
(I2.16P is necessary and sufficient for stability of the discretization of the Stokes equation 
with homogeneous boundary conditions. The inf-sup condition (12.171) and the restriction 
(I2.18P are standard conditions for stability of the discretization of the Darcy equations 
with homogeneous boundary condition on the normal trace. 

Condition (I2.20p is the existence of a uniformly bounded right-inverse of the operator 
Ho(f2o) 3 V/j Vft • n G $5?. A s discussed in [2U Section 5], this condition is satisfied for 
Brezzi-Douglas-Marini (BDM) and Raviart-Thomas (RT) elements (see below for their 
definitions) on general shape-regular triangulations in the plane, and on tetrahedrizations 
of the space that are quasi- uniform near the boundary S. Existence of satisfying (12.201) 
for BDM and RT elements in general tetrahedrizations is an open question. Hypothesis 
(I2.19P is a very mild condition demanding that the discrete space for the Stokes condition 
can provide non-trivial flow in to the Darcy domain without a blow-up of the velocity 
field. This condition is discussed in [241 Section 6], where it is shown that as long as the 
Stokes velocity space contains piecewise linear functions, this condition is satisfied. 
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Some examples For precise descriptions of the finite element spaces below, the reader is 
referred to [7J, [17] and [19] . All choices below will be given with the following assumptions: 

• Hypothesis (12.191) will be assumed to hold. 

• The convergence orders of the Stokes and Darcy elements are chosen to match. 

• If the discretization is conforming, we will assume that the Darcy partition E^ is 
either equal to or a refinement of Eg. 

The Brezzi-Douglas-Marini (sometimes called Brezzi-Douglas-Duran-Fortin in the three 
dimensional case) is the mixed finite element that uses the spaces 



for k > 1. We will refer to it as the BDM(fc) element. The BDM(l) element can be 
coupled in a conforming way with the MINI element and the Bernardi-Raugel element. It 
can also be coupled with the P2-iso-Pi element in a conforming way if E^ is either equal 
to or a refinement of E s and in a non-conforming way otherwise. The BDM(2) element 
can be coupled in a conforming way with the conforming Crouzeix-Raviart element. More 
generally speaking, BDM(fc) can be coupled with the Taylor-Hood element of order k for 
any k > 2. 

The Raviart-Thomas element of order k, henceforth referred to as RT(fc), is defined 



where RT fc (T) = {p(x) + g(x) x : p G F k (T) d , q G P fc (T)}. The RT(0) element can be 
coupled in non-conforming way with the MINI element and the Bernardi-Raugel element. 
For k > 1, RT(/c — 1) can be coupled with the Taylor-Hood element of order k. 

3 An alternative point of view 

In this section we propose a different way of interpreting the coupled method, based on 
seeing the Darcy equations as part of a generalized boundary condition for the Stokes 
problem. 

3.1 The Darcy boundary condition 

Given /d satisfying the compatibility condition (I2.5p . we consider the solution of the 
Darcy problem 



H h (fi D ) := {u fe GH(div,fi D ) : u h \ T e¥ k (T) d VT G 7#}, 
L h (Q D ) := {p h :n D ^R : Ph \ T eF k ^(T) VT G T D h }, 



as the pair 



H h (Q B ) := {u h GH(div,fi D ) : u h \ T eRT k (T) VT G T£} 
L h (n D ) := {p h : Q D M : p h \ T G P fc (T) VT G T D "}, 



K + = in Q D , 
divu£ = / D in Q D , 
Up • n = on Td U E, 
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and note that p^ G if^^o)- Also, for <p G £ 2 (£), we consider the solution of 



■D + V % 


= 


in f^. 






in Qd. 




/* 

It, 


U D n 


= 


on r D 


<t> 

U D n 


= 


on E, 


/ vi 


= o, 











and define with it the Flux-to-Pressure operator FtP(0) := p^ls- It is simple to prove that 
FtP is a linear and symmetric operator in L 2 (S). Indeed, the Flux-to-Pressure operator 
satisfies 

(FtP(0), v • n) s = a D « v) - 6 D (v,p D ) (3.1) 
for all v G H D (div, f2 D ) such that v ■ n G £ 2 (£), which gives 

(FtP(0),V>£ = (FtP(0),uS -n) s = a D (u£,ug) = (FtP(^),0) 2 V0,^ G L 2 (E). (3.2) 
It is clear that, as L u s • n = 0, we have the splitting 

P D |s=^ + FtP(u s -n) 

for the Darcy pressure on S. This allow us to write the coupling conditions (12. 4 j) as a 
unilateral boundary condition for the Stokes flow on the interface E: 

2z/e(u s )n - p s n + K7r t u s + FtP(u s ■ n)n = -pLn. (3.3) 
v v ' 

The underbracketed term corresponds to a symmetric positive semidefinite non-local op- 
erator that takes into account the influence of the Darcy domain on the Stokes flow, acting 
separately on the tangential and normal components of the Stokes flow. The Stokes sys- 
tem (12 .ip can then be complemented with the non-local condition (13. 3 j) and the Dirichlet 
condition on Ts (see (12. 3p ) to produce a formulation of the Stokes-Darcy problem that is 
equivalent to (I2.13p . It consists in looking for ug G Hg(f2g) and ps G L 2 (f2g) such that 

a s (us,v) + c(u s ,v)-6 s (v,p s ) = £(v) Vv G H^(O s ), 

(3.4) 

6 s (u s ,g) =0 VgGL 2 (fi s ), 

where 

c(u,v) := (FtP(u-n),vn) s (3.5) 

and 

£(v) := (f s ,v) Qs - (p{>,v -n) s . 

By (13.21) . it follows that the bilinear form in (13.51) is symmetric and positive semidefinite. 
A simple argument shows that the bilinear form as(us, v) + c(ug, v) is coercive in H 1 (fis)- 
This fact gives a very simple proof of the fact that the Stokes-Darcy system is well posed 
and that it can be understood as a modified Stokes problem without losing any of its good 
properties. This will be exploited to design an effective Krylov-based iterative method 
to solve the algebraic linear system of equations arising from the discrete counterpart of 

(E3D. 
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3.2 The discrete flux-to-pressure operator 

If we now choose discrete spaces for the Darcy problem satisfying (12.17p - fl2.18j) . we can 
define a discrete version of the operator FtP as follows. Given 0^ G $^ with J E <ph — 0, 
we define FtP/ l (0/ l ) : $p — > M. to be the functional (compare with (13.1 1) ) 

(FtP ft (0 h ), v h • n) E := a D « v h ) - & D K,rf), Vv„ G H*(fi D ), (3.6) 
where (u^,p^) G H^^d) X Lq(^d) solves the discrete equations: 

u h ■ n = 4> h on £, (3.7) 

u ft • n = on T D , (3.8) 

a D (u£v fc ) -bn(v h ,pt) = Vvft G H£(fi D ), (3.9) 

6o(ut?fc) = Vg, GiJ(0 D ). (3.10) 
With arguments similar to those used in the continuous case, it is easy to prove that 
(FtP ft (0 A ), t/Jhh = a D (uf , uj) = (FtPftC^h), <Ms V^, ^ G 

which shows that the discrete flux-to-pressure operator FtP^ is also symmetric and non- 
negative. 

The discrete pressure due to sources, j[, can be similarly defined as a residual: 

where (u{,p{) G Hq^d) x -^o(^d) solve the discrete equations: 

a D (u£,v A ) -b D (v h ,p{) = Vv^gH^(Od), 

br>(u f h ,q h ) = (/ D ,gfc)oD Vg h , G I#(fi D ). 

We recall that the operator is the L 2 (S)-projection on It is straightforward 
that the discrete Darcy pressure and velocity of problem (12.151) admit the splitting 

Pn=Ph+Ph and u£ = u£ + V Vs ; . 

It follows that (12.151) may be equivalently stated as follows: find (ug,p|) G Vg(fig) x 
L fo (fi s ) such that 

a s «v ft ) + c ft (u§,v fc )-6s(v A ,^) =4N Vv h eVt(n a ), 

(3.11) 

&s«g h ) =0 Vg fc e L fc (fis), 

where 

c h (u h , v A ) := (FtP ft (i2 ft (u h • n)), i^(v h ■ n)) E 

and 

4K) := (f s , v h )n s - (P f h , R h (v h • n)) E . 

Inn the conforming case ($g C the L 2 (S)-projection operator R h does not play any 
role in the formulation. 
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4 The decoupled iterative method 

We introduce the self-adjoint operators A$ and Cg denned from Vg(fig) to its dual 
Vg(fis)*by 

(A s u, v) V h (Qs) , xV h (ns) =a s (u,v) and (CgU, v) V |(n s )*xV|(n s ) = c /i( u > v ) 
for all u, v G Vg(O s ). Let also Bg 1 : V§(ft s ) ->- L h (fi s )* be the operator defined by 

(BsU, q) L h (n s )*xL h (n s ) = _ ^s(u, q) 

for all u G V§(ft s ) and g G L fe (fi s ). 

Problem (13. lip can be written in operator form as follows: 

where 

A * '■= if*B^ {Bk f) '■ V ^ s) X Lk{ns) V ^ s) * X 

and (-Bg)* is the adjoint of 5g. We know from Theorem 12. II that both ||-4g|| and ||(^4g) _1 || 
are uniformly bounded in h. If we denote by /g : L h (Qs) — > L h (Qs)* the Riesz operator 
defined by 

(!sP> q)Lh(n s )*xLh(n s ) = (P, q)n s Vp, q G L h (Q s ), 
then, the positive-definite self-adjoint operator 

rl ;= (f ^ : v*(n s ) x l*(o s ) -> v h s (n s y x L h (n s y 

and its inverse are uniformly bounded uniformly in h. It follows that the condition number 
of (Pg) _1 ^lg is bounded from above by a constant independent of the mesh parameter 
h. Consequently, the MINRES algorithm preconditioned with (Pg) -1 solves ( 14. ip with a 
reduction of the norm of the residual that is independent of the mesh size h. 

Let us now discuss how action of Cg on a given Ug G Vg(fig). To this end we introduce 
the self-adjoint operators A^ and defined from Hq^d) to its dual Hq(^2d)* by 

(^d u . v )h^si d ).xh^(si d ) = Od(u,v) and (^u,v) H k (!!D) , xHS(S j D) = (divu, div v) nD 

for all u, v G Hq(Q d ). Let us also consider B^ : Hq(Q d ) -> Lq(Q d )* given by 

(B D u, g)^(n D )*xi§(n D ) = - & d(u, q) 

for all u G H£(fi D ) and q G Z#(fl D ). 

We compute Cgiig through ( 13.6f) after solving problem ( 13. 7ft with 0^ = ii^Ug ■ n). 
This is to say that we have to deal with a saddle point problem of the form 

A (8 = S) (4 - 2) 

n 



where 



and (B^Y is the adjoint of Bfc. 

The stability of the pair of spaces (H{f(fi D ), Lg(O s )) fl2TT7D - fl2TT8j) ensures that both 
||.AdI| and 1 1 ) 1 1 1 are uniformly bounded in h. If we denote by 1^ : Lg(fi D ) — >• Lq(Q d )* 
the Riesz operator given by 



( I DP^)L' i (n D )*xL^(n D ) = (p,q)n VP, g e L (fi D ), 
it is clear that the block diagonal positive-definite self-adjoint operator 

n := + D » ^ : H*(n D ) x L h (n D ) -> H^ D )* x L^ D )* 

and its inverse are also uniformly bounded in ft,. It follows that we can find an inclusion set 
for the eigenvalues of ("Pd) _1 Ad that is independent of h. This means that the MINRES 
method preconditioned with (V^)~ l yields the solution of (14.21) in a number of iterations 
independent on the mesh size h. 

Summing up, the decoupled iterative method we are proposing here to solve ( I2.15P 
consists in two nested MINRES algorithms. Computationally, the actions of the precon- 
ditioned correspond to solving two decoupled local problems. The first one is defined 
by the bilinear form as(v) m Vg(fig) and corresponds to the block Ag. Actually, A§ 
is associated with the operator — 2i/div(e(-)). Therefore, the local problem in the fluid 
amounts to a vector Laplace equation with a Dirichlet boundary condition on T$, a Neu- 
mann condition in the normal direction and the slip boundary condition in the tangential 
direction on S. The other local problem is defined by the bilinear form 

(K _1 u D ,v D )n D + (divu D ,divv D )n D 

on Hq(^ d ) corresponding to the diagonal block A^ + D^. 

For the construction of practical preconditioners for discrete systems, the computa- 
tional cost of evaluating these operators and the memory requirements of these procedures 
are key factors. The exact inverses appearing in the canonical preconditioners should be 
replaced by proper cost effective, and norm equivalent operators. Let us consider self- 
adjoint and positive-definite operators Pg and P£ that are spectrally equivalent to A^ 
and A^ + respectively, i.e., 

(AgU S , U S ) Hs (Q s )*xH s (f7 s ) - (-Ps U S, u s)h s (^s)*xH s (Q s ) 

and 

((A D + -D D )u D , U D ) Ho (fi D )*xHo(n D ) - ( P D U D, Ud)h (S] d )«xH (S] d )- 

Then, instead of {V^)~ l and we can use respectively the preconditioners 

0P s V o \ , o 



o (il)- 1 ) and V o 
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and still have an optimal decoupled iterative method for problem f)2.15p . Ideally, we 
would have the actions of (Ps - ) -1 and (Pfj) -1 cos ^ about the same as the actions of Ag 
and + D^. As A$ corresponds to a second-order elliptic operator in Hg(fig), we can 
easily take advantage of multigrid techniques or domain decomposition methods to find a 
good candidate for (Pg 1 ) -1 . The construction of a preconditioner (P^) -1 is l ess obvious. 

4.1 Nodal auxiliary space preconditioning in H(div) 

In this section we describe the construction of the nodal auxiliary space preconditioning 
of Hiptmair and Xu [20] for elliptic problems in H (div, Q-d). This is our choice here for 
the matrix version (P^) -1 of the preconditioner (P^) -1 needed in the last section. To fix 
the ideas, we assume that the tensor K is given by t~ 1 \ where I is the identity in M. d and 
t is a given positive constant. In our numerical experiments, Hq^d) is derived from the 
KT(k — 1) or BDM(fc) mixed finite elements with k = 1 or 2. Let i = 1, ...,/} be 

the usual basis of Hq^d), then the matrix realizations of A 1 ^ and are given by 

A D : = { T {<t>i,<i>i)nnh<i,j<I 

and 

d d : = ((div^div^jnji^jx/ 

respectively. Our aim is to provide a matrix 6l M that is spectrally equivalent to 
+ Dp and such that the action (P^) -1 on a given vector is easier to compute then 

that of (A£ + D^)- 1 . 

We denote by [V h (n D )] d C Hj(fi D ) the standard space of piecewise and continuous 

vector fields and consider its usual nodal basis {<p e , £ = 1, • • • ,Ld}, where L is the 

dimension of V^^d)- We introduce the matrix given by 

(L ft )4* = {V<p e , V<p m )n D + r{(p £ , <p m )n D , l<£,m<Ld. 

In the three dimensional case (d = 3), we also need to consider the Nedelec space 
WJ(Od) C H (curl, f2o) of order k. We denote its usual basis {<Tj, i — 1, . . . ,N}. We 
introduce the diagonal matrix S^j url given by 

(Sr% := (curl en, curl i — 1, ■ ■ ■ , N 

and denote the diagonal of A^ + by S^ lv . 

In the three-dimensional case, we denote by G M iYx/ the matrix that represents 
curl : Wj(n D ) — y Ho(fi D ) in the following sense, 

i 

curler, = y^Cftjjj^j-, Vz = 1, • • • , N. 

In the the two-dimensional case, the matrix G M Lx/ is defined similarly with respect 
to the operator curl : Vq(Q-d) — > Hg(f2p) defined by curlw = 
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We use II™ and II^ 1V to denote the canonical interpolation operators onto the finite 
element spaces Wj(0 D ) and Ho(f2 D ) respectively. The mappings Il^ lv : [V^^d)^ — > 
Hq(Qd) and n™ rl : [V^d)] 3 ->■ W£(fi D ) will be described by the matrices lf v G R dLxI 
(d = 2, 3) and I£ url G R 3LxN defined by 

i 

3=1 

and 

respectively. 

The 3d-H(div) auxiliary space preconditioner of Hiptmair and Xu consists in 

(p^- 1 := (sf j- 1 +if(L ft )- i (ir) t +r- i c ft ((sr 1 )- 1 +ir 1 (L fc )" 1 ar 1 ) t )cji 

and the 2-d version of this preconditioner is defined by 

(PdT 1 := (Sf)- 1 + lf(L /l )- 1 (ir) t + r^C, (-A,)" 1 C* 

where the matrix A^ stands for the discrete Laplacian on the finite element space Vq(Qd). 

Notice that the transfer matrices C^, I^ 1V and I^ url corresponding to the curl operator 
and the interpolations are sparse matrices that can be computed in a straightforward 
manner. The evaluation the preconditioner is then essentially reduced to several second- 
order elliptic operators. Hence, standard multigrid techniques domain decomposition 
methods for H l equations can be applied. 

5 Numerical results 

This section is devoted to the description of numerical experiments validating the effec- 
tiveness of the decoupled iterative method. We will show results for two dimensional 
problems, considering three examples of pairs of stable elements for the Darcy-Stokes 
problem. The first two examples correspond to the conforming Galerkin schemes based 
on the combination of the MINI and P 2 -iso-Pi elements for the Stokes problem with the 
lowest order Brezzi-Douglas-Marini element BDM(l). The third one is the nonconform- 
ing scheme resulting from the Taylor-Hood element and the second order Raviart-Thomas 
element RT(1). 

5.1 Convergence rates 

We begin by introducing some notation. The variable DOF stands for the total number 
of degrees of freedom defining the finite element subspaces X h and Q h , and the individual 
errors are denoted by: 

e(u D ) := ||u D - UD|| H (div,f} D )> e ( u s) := ll u s - Ug|| H i(fi s ), 
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and 



where 



e(PD; 



\Pd 



rr) 1 II 



e(ps) : = 
2d and pg 



IPs 



Pslln s , 



u fe |n D , Ug := u/i | si D , := P^d ancl P's := Ph\si s witn ( u ^pO G X x Q 
being the solution of (|2.15|) . We also let r(u D ), r(u s ), r(p D ) and r(p s ) be the experimental 
rates of convergence given by 



and 



r l u D, 



t(pd) 



log(e(u D )/e / (u D )) 
\og{h/h') 

log(e(p D )/e'(pD)) 



r l u sJ 



r(ps 



log(e(u s )/e'(u s )) 
log(h/h') 

\og(e(p s )/e'(p s )) 
\og{h/hf) 



\og(h/h>) 

where h and h' are two consecutive mesh sizes with errors e and e'. 

We now describe the data of the example. We consider the domains Qd '■— (0, 1) x 
(0, 1/2) and fi s := (0, 1) x (1/2, 1). We take v = 1, k = 1 and K = I, the identity of 
IR 2x2 . The right hand side functions are selected in the model in such a way that the 
exact solution is given by: 



Pd(x) 



6tt 



%2 1 

sin(27ra;2 

2 4vr v 



sin 2 (27nri) cos(27ra;i), 




sin(7rx 2 ) cos(7re 2 ) sin 3 (27ra;i) 
-3 sin 2 (27rxi) cos(27rxi) sin 2 (7rx 2 ) 



ps{x) : = cos (7^1 J ^2 + 0.5-2 cos 2 (7^2 + 0.5) 
in the Stokes domain. 



DOF 


h 


e(u s ) 


r ( u s) 


e(ps) 


r(ps) 


e(u D ) 


r(u D ) 


e(PD) 


r (PD) 


543 


1/8 


1.86E+01 




9.26E-00 




4.73E+01 




1.60E-01 




2043 


1/16 


1.01E+01 


0.87 


3.04E-00 


1.60 


2.48E+01 


0.92 


8.10E-02 


0.98 


7923 


1/32 


5.17E-00 


0.97 


8.80E-01 


1.79 


1.26E+01 


0.98 


3.99E-02 


1.02 


31203 


1/64 


2.59E-00 


0.99 


2.49E-01 


1.82 


6.31E-00 


0.99 


1.98E-02 


1.00 


123843 


1/128 


1.29E-00 


0.99 


7.56E-02 


1.72 


3.16E-00 


0.99 


9.92E-03 


1.00 



Table 1: Convergence rates: MINI-BDM(l) 



We begin by providing a numerical exploration of the asymptotic convergence rates 
of the three examples. In Tables [Q [2] and El we summarize the convergence history of 
the Galerkin scheme (I2.15P for a sequence of nested uniform meshes of the computational 
domain Q := (0, l) 2 by means of triangles. All the results are obtained by applying our 
decoupled preconditioning technique. In each case we display the numerical rates of con- 
vergence versus the degrees of freedom DOF. We observe that, as expected, in the case 
of the MINI-BDM(l) and the P2-iso-Pi-BDM(l) couplings, the convergence is linear for 
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DOF 


ft 


e(u s ) 


r ( u s) 


e(ps) 


r(ps) 


e(u D ) 


r(u D ) 


e(PD) 


r (PD) 


385 


1/8 


1.86E+01 




4.10E-00 




4.73E+01 




1.60E-01 




1423 


1/16 


1.01E+01 


0.87 


2.14E-00 


0.94 


2.48E+01 


0.92 


O 1 1L 1 AO 

8. HE— 02 


n no 

0.98 


5467 


1/32 


5.17E-00 


0.97 


6.03E-01 


1.82 


1.26E+01 


0.98 


3.99E-02 


1.02 


91 /L97 


1 /P.A 


2.59E-00 


0.99 


1.57E-01 


1.94 


6.32E-00 


0.99 


1 QSP 09 


l.UU 


O/IQOC 

o4oo0 


1 / 12© 


1.30E-00 


0.99 


4.25E-02 


1.88 


3.16E-00 


0.99 


y.yzhv— Uo 


l.UU 






Table 2: Convergence rates: P2 


-iso-Pl-BDM(l) 






DOF 


h 


e(u s ) 


r ( u s) 


e(ps) 


r(ps) 


e(u D ) 


r(u D ) 


e(PD) 




887 


1/8 


4.09E-00 




1.08E-00 




1.48E+01 




5.23E-02 




3371 


1/16 


9.56E-01 


2.09 


8.88E-02 


3.60 


4.03E-00 


1.87 


1.35E-02 


1.95 


13139 


1/32 


2.37E-01 


2.01 


7.06E-03 


3.65 


1.07E-00 


1.90 


3.40E-03 


1.99 


51875 


1/64 


5.93E-02 


1.99 


7.85E-04 


3.17 


2.94E-01 


1.87 


8.50E-04 


2.00 


206147 


1/128 


1.48E-02 


1.99 


1.98E-04 


1.98 


8.43E-02 


1.80 


2.12E-04 


2.00 



Table 3: Convergence rates: Taylor-Hood-RT(l) 



the velocities in both the Stokes and the Darcy domains. The Taylor-Hood-RT(l) scheme 
provides a quadratic convergence for the Stokes and Darcy velocity unknowns. We fixed 
the tolerance parameter for the outer MINRES method at 10~ 6 and checked empirically 
that the largest inner MINRES tolerance parameter that provides a convergence in agree- 
ment with the rates predicted by the theory is 1CT 2 . All the results displayed here are 
obtained with this combination of tolerance parameters. 

5.2 Performance of the iterative method 

In the following, we will denote by Ag, Bg, Cg and Mg the matrix realizations of A$, 
.Bg, Cg, and 7g respectively. Similarly, A^, B^, and are the matrix realizations 
of Bp, Dp, and 1^ respectively. 

The numerical results were obtained using Matlab's own MINRES routine. For all 
experiments, the convergence is attained when the Euclidean norm of the relative residual 
is reduced by 10~ 6 for the outer MINRES while (as indicated above) the tolerance for the 
inner MINRES method is set to 10~ 2 . The outer MINRES is applied to a linear system 
of equations with matrix 

/Ag + C| (Bgn 

v B s o 

It is initialized with the solution of the Stokes problem with a non slip boundary condition 
Tg and an homogeneous Neumann boundary condition on S. The MINRES algorithm is 
accelerated with one of the following preconditioners: 

/(Ag); 1 \ BPX ._ /(A^-p 1 , \ 

/s -~V o (Mi)- 1 ;' /s ~ v ( M s)~V' 

where the notation (Ag); 1 means that the linear systems of equations with matrix Ag are 
solved by a direct solver (with the Matlab backslash command) while (Ag)^ means that 
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we use the Bramble-Pasciak-Xu [6j [28] preconditioner corresponding to the SPD vectorial 
Laplace matrix Ag. 

On the other hand, each application of Cg to a vector requires the solution of a saddle 
point problem with matrix 

(B£>) tN 





We again accomplish this task applying the MINRES method preconditioned with one of 
the following symmetric and block diagonal matrices: 



;a& + d&) 









v; 











BPX 



,Pd)bpx 








In the definition of the preconditioner V^, (A^ + Dp)" 1 means that we simply use a direct 
solver for A^ + with the aid of the backslash Matlab command. The preconditioners 

and 7-0 PX are obtained by substituting (A^ + D^) -1 in by the Hiptmair and Xu 
preconditioner (Pj?,) -1 . The subscript \ in (P^,)" 1 means that we solve the underlying 
Laplace problems with a direct solver, with the Matlab backslash command, and (P^bpx 
means that we use the well-known BPX-preconditioner [6], [28] for (Lft) -1 and (—Aft) -1 . 

In the cases where the mass matrix is diagonal the action of its inverse can be explicitly 
computed. In the other cases, one simple and effective strategy consists in substituting 
the action of the inverse of the mass matrix by one sweep of the symmetric Gauss-Seidel 
method. 



DOF 


h 












^BPX^BPX) 


543 


1/8 


26(4) 


26(26) 


26(29) 


56(4) 


56(26) 


56(29) 


2043 


1/16 


32(4) 


32(30) 


32(46) 


84(4) 


84(30) 


84(46) 


7923 


1/32 


40(4) 


40(33) 


40(62) 


121(4) 


121(33) 


121(62) 


31203 


1/64 


46(4) 


46(38) 


46(77) 


144(4) 


144(38) 


144(77) 


123843 


1/128 


50(4) 


50(42) 


50(91) 


158(4) 


158(42) 


158(91) 



Table 4: Number of iterations: MINI-BDM(l) 



DOF 


h 












^BPX^BPX) 


385 


1/8 


24(4) 


24(26) 


24(29) 


50(4) 


50(26) 


50(29) 


1423 


1/16 


30(4) 


30(29) 


30(46) 


80(4) 


80(29) 


80(46) 


5467 


1/32 


36(4) 


36(34) 


36(62) 


107(4) 


107(34) 


107(62) 


21427 


1/64 


42(4) 


42(38) 


42(76) 


130(4) 


130(38) 


130(76) 


84835 


1/128 


44(4) 


44(41) 


44(91) 


146(4) 


146(41) 


146(91) 



Table 5: Number of iterations: P2-iso-Pl-BDM(l^ 



In tables HJ [5] and EJ we list the number of iterations of the two nested MINRES 
methods with different combinations of preconditioners. The preconditioner in brackets 
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DOF 


h 






887 


1/8 


28(5) 


28(28) 


3371 


1/16 


34(5) 


34(32) 


13139 


1/32 


38(5) 


38(36) 


31203 


1/64 


42(5) 


42(40) 


206147 


1/128 


44(5) 


44(44) 



Table 6: Number of iterations: Taylor- Hood-RT(l) 



is the one used for the inner MINRES. We show the number of outer MINRES iterations 
and the number in brackets is an average of the number of inner MINRES iterations. We 
observe that for different mesh sizes, the iterative method results in a uniform number of 
MINRES iterations. Therefore, the preconditioners are robust with respect to the mesh 
size, which agrees with the theoretical results. 
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